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We report on a non-perturbative approach to the ID and 2D Hubbard models that is capable of 
recovering both strong (U S> t) and weak-coupling (U <C t) limits, with U the on-site Coulomb repulsion 
and t the kinetic energy. Dynamical corrections to the electron self-energy in the single particle Green 
function are explicitly included by expanding in terms of the 16 eigenstates that characterise two 
nearest neighbour sites. We first show that even when U is much smaller than the bandwith, the 
Mott-Hubbard gap never closes at half-filling in both ID and 2D. Consequently, the Hubbard model at 
half-filling is always in the strong-coupling non-perturbative regime. For both large and small U, we 
find that the population of nearest-neighbour singlet states approaches a value of order unity as T — » 
as would be expected for antiferromagnetic order. We also find that the double occupancy is a smooth 
monotonic function of U and approaches the anticipated non-interacting limit of 1/4 as U — > and 
vanishes as U — > oo. Finally, we compute the heat capacity (C(T, U)) for both ID and 2D. Our results 
for ID at moderate to high temperatures are in quantitative agreement with those of the exact Bethe 
ansatz solution, differing by no more than 1%. In addition, we find that in 2D, the C(T,U) curves 
vs T for different values of U exhibit a universal crossing point at two characteristic temperatures, 
T ~ 1.7t ± 0. It and T ~ 0.4 ± O.lt as is seen universally in Hubbard models and experimentally in a 
wide range of strongly-correlated systems such as 3 He, UBes, and CeCu%- x Al x . The success of this 
method in recovering well-established results that stem fundamentally from the Coulomb interaction 
suggests that local dynamics are at the heart of the physics of strongly correlated systems. 



I. INTRODUCTION 

In both the weak and strong coupling regimes, the 
Hubbard model at half-filling is expected to be an an- 
tiferromagnet at T = 0. The strong-coupling argument 
relies on an isomorphism between the half-filled Hubbard 
model in the limit that the on-site Coulomb repulsion, U, 
exceeds the hopping integral, t, and a Heisenberg antifer- 
romagnet. As a consequence, the corresponding ground 
state is antiferromagnetic and the energy scale for the 
gap between the upper and lower Hubbard bands is set 
by the energy cost for double occupancy, A « U. In 
the opposite or weak-coupling limit, U <C t, perturba- 
tion theory predicts that a van-Hove singularity induces 
a spin-density wave producing a gap that is exponentially 
small in the Coulomb repulsion; that is, A « ie -27I V t / !7 . 
While the perturbative argument cannot be extended 
into the strong-coupling regime where the t/U mapping 
to a perfect Heisenberg antiferromagnet applies, conti- 
nuity between the two regimes suggests that antiferro- 
magnetism persists for any non-zero value of U. Further, 
this argument would also suggest that the Hubbard gap 
never closes for any non-zero value of U. Consequently, 
the half-filled Hubbard model is always in the strong- 
coupling regime. However, no exact results are known. 
In fact, while numerical simulations [Q support an anti- 
ferromagnet at half-filling, several mean-field arguments 
suggest otherwise. Dating back to the pioneering work of 



Mott |^] and Brinkman and Rice numerous calcula- 
tions on the 2D @,|] or the D = oo ||| half-filled Hub- 
bard models suggest that whenever U is much smaller 
than the bandwidth, W — 8t, the Hubbard gap closes 
and a metallic phase ensues. Contrastly, Anderson [|| has 
argued that almost certainly the half-filled 2D Hubbard 
model is non-perturbative as in the ID case ||, thereby 
possessing a discontinuity only at U = 0. At the heart 
of the non-perturbative nature of the Hubbard model g] 
is the projective mismatch between the low-energy phys- 
ical subspace and the "anti-bound states" which form in 
2D for any non-zero value of U. Antiferromagnetism fol- 
lows necessarily as a corollary from the break-down of 
perturbation theory. 

In this paper we re-examine this problem using an ap- 
proach that is capable of spanning the weak and strong- 
coupling regimes. Our approach is based on the Hub- 
bard operators which exactly diagonalize the interac- 
tion part of the Hubbard Hamiltonian. Consequently, 
the Hubbard operators are tailor-made to access the 
strong-coupling regime, U 3> t. Rather than work in 
the static approximation in which quantum fluctuations 
are ignored leading to infinitely sharp upper and lower 
Hubbard bands, we include the dynamical corrections 
which lead to broadening of the spectral features. As 
in the work of Matsumoto and Mancini Jlt^U] , we fo- 
cus on the dynamics associated with two neighbouring 
sites. The dominant dynamics appear to be governed 
by spin-fluctuations which lead to singlet-triplet excita- 
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tions. From our analysis, we conclude that the Hubbard 
gap never closes and the 2D Hubbard model at half-filling 
is always in the strong-coupling regime. Our results then 
corroborate those of a recent improvement fL2| on dy- 
namical mean-field theory |7j] in which the momentum- 
dependence of the self-energy is explicitly included |E| at 
particular points in the Brillouin zone. To determine the 
validity of the approach we use here, we study as well the 
ID half-filled Hubbard model as exact results are known 
from Bethe ansatz [JDS. As expected, we find that the 
Hubbard gap persists even in the weak coupling regime. 
In addition, we find that our results for the heat capac- 
ity are in perfect agreement with those from Bethe ansatz 
in the temperature range where the Coulomb interaction 
dominates the physics, that is, moderate to high tem- 
peratures. Finally, we show that we recover the well- 
established universal crossing of C(T, U) vs T for 
various values of U that is seen experimentally in a wide 
range of strongly-correlated systems such as 3 He [17|, 
CeCu 6 _ x Al x @, Ndi- x Ce x Cu0 4 @, and UBe 3 @j. 
The success of our approach suggests that local dynam- 
ics lead to many of the features of strongly-correlated 
electronic systems. 



II. DYNAMICAL GREEN FUNCTION 
APPROACH 

The starting point of our analysis is the on-site Hub- 
bard model 



H 



tin Ca —C-i 



l i\ n il 



(1) 



where Uj = t if are nearest-neighbour sites and 

zero otherwise. Rather than working with the origi- 
nal electron operators, we use the Hubbard operators 
Via = c ta ni^ a and = c icr (l - nj_ CT ) as these opera- 
tors exactly diagonalise the interaction term. In terms of 
the Hubbard operators, Ci„ = Via + Cia- While the inter- 
action term is now simplified in this basis, the Hubbard 
operators do not obey standard Fermi statistics, making 
impossible any diagrammatic approach based on Wick's 
theorem. However, the equation of motion approach has 
been demonstrated ||l]-|23| to offer an alternative to the 
diagrammatic expansion. Consider the two-component 
basis 



Si<7 
Via 



(2) 



and its associated Green function G(i,j,t,t') = 

((V^;V4» = ^-«0({^(i),f] CT (t')}) where {A, B} 
is the anticommutator and (• • •) is the thermal average. 
The equations of motion for the Hubbard operators, 



.OA 
'-dt 



(3) 



will of course contain a contribution which is linear in the 
Hubbard basis and in addition new terms, Sj® which con- 
tain operators that lie outside the Hubbard basis. Ideally, 
if such operators are included in the Hubbard basis, then 
the non-linear contributions can be minimized. However, 
such a procedure is necessarlily cumbersome. Instead of 
enlarging the basis, we project 8jf onto the Hubbard ba- 
sis using the Roth p3| projector, 



V{0)=Y d {{OM))Kn^n 



(4) 



In 



which projects any operator O onto the Hubbard opera- 
tor basis, where 7(k) = FT{{il>„{i, t), ipl(j, t}), and FT 
denotes the time and space Fourier transform. This pro- 
jector is particularly useful because it allows us to recast 
the equations of motion for the Hubbard operators 

ji(t) = E A + V(6j?) + 8ji = + Sj z (5) 

in terms of the renormalized energy matrix 

E(k) =E + FT({5jl^i})I- 1 (k) (6) 

and a correction 5ji = 5j° — V(Sjf). Clearly, V(5ji) = 
0. The formal solution for the Fourier transform of the 
Green function 



G(k,w) 
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uj- E(k)- Sm(k,uj) 
contains the dynamical self-energy, 



I(k) 



5m{k,uj) = FT{R{5j ia {t),5jl{t')}) I 



(7) 



(8) 



where R denotes retarded and I the irreducible part. For 
a paramagnetic phase, the overlap matrix 



1-| 

o i 



h o 

h 



(9) 



is explicitly diagonal. The weights which appear along 
the diagonal represent the contribution from the lower 
and upper Hubbard bands, respectively. Note, they sum 
to unity. This feature coupled with the fact that the dy- 
namical corrections vanish when U — guarantees that 
we recover the correct non- interacting limit. 

The primary operational hurdle with any analytical ap- 
proach to the Hubbard model is the evaluation of the dy- 
namical self energy. In the static approximation [pl| p3[, 
the self-energy is dropped, and the Green function re- 
duces to the pole structure 



G(k,w) 



-E(k) 



J(k) 



(10) 



where -E'(k) defines the energy bands. At this level of 
theory, the Hubbard bands are sharp as the Green func- 
tion has a pole at the energy of each band. As our focus, 
however, is on the closing of the Hubbard gap, the pole 
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approximation is inadequate and the broadening in the 
bands arising from the self energy is crucial. In the con- 
text of the composite operator approach, Mancini and 
Matsumoto JlCj have developed a real-space scheme for 
computing the dynamical corrections to the static ap- 
proximation. To implement this procedure, we rewrite 
the dynamical correction 



<5m(k, ui) = Dm()i,Lu) 



1 -1 
-1 1 



(11) 



FT(RSji(t)6jJ(t))K 



in terms of the 2x2 matrix K. Because Dm(k, oj) cannot 
be evaluated exactly, we seek a systematic way of calcu- 
lating the dynamical corrections. The simplest approach 
would be to consider the single-site approximation. Such 
an approximation is in the spirit to the d = oo || meth- 
ods in which the self-energy is momentum independent. 
An improvement would be to consider the dynamics as- 
sociated with two sites as proposed by Mancini and Mat- 
sumoto [ ^Opl| ] . Evaluation of the self-energy over succe- 
sively larger clusters would lead to an exact determina- 
tion of the dynamical corrections. Hence, we write the 
dynamical corrections as a series 

Dm(x, x') = 8 XjX iDmo{x,x') + y] S x+atX >Dmi (x, x') + ■ ■ 



in increasing cluster size. Here, x and x 1 are neighbour- 
ing sites and a indexes all nearest-neighbour sites. In 
the two-site approximation, the series is truncated at the 
level of on-site, Dmo, and nearest-neighbour, Dmi con- 
tributions. In Fourier space, the dynamical corrections 
can be written as, 



with 



Dm(k, lu) Dmo(uj) + a{\f.)Dmi{ui) 



Dm Q (Lu) = -FT(R{5j(t),5]Ht')}), 
Dm 1 (u J ) = ^FT(R{Sj(t),Sf\t')}) 1 



(12) 
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where the factor of 1 /4 arises from the coordination num- 
ber on a square lattice, a(k) = (cos k x -t-cos k y )/2, and Sj 
and Sj' are centered on x and x' , respectively. The goal 
here is to study the half-filled Hubbard model by evalu- 
ating Dmo and Dm\ in the two-site approximation. 

Because the procedure for implementing the two-site 
approximation has been described previously ]Io| , [rH , we 
will only outline the essence of the method: 1) enumer- 
ate the quantum mechanical states for two neighbouring 
sites, 2) use the resolvent method [|4| to express how the 
surrounding environment interacts with the two-site sys- 
tem, 3) express the Sj's in terms of the level operators 
for the two-site system, and 4) use the non-crossing ap- 
proximation to expand Dm(k, u>) in terms of the two-site 



states. For two sites, there are 4 2 = 16 electronic states. 
Let B(i), F a (i), and D(i), represent single site level oper- 
ators acting on empty (a boson state), singly (with spin 
a) and doubly occupied sites respectively. In terms of 
the level operators, the original Hubbard operators are 
£er = B^Fc and r\ a = aF_ a D. The level operators for 
two-site states, <!>„, are formed by all possible symmetric 
and antisymmetric combinations of B, F a and D on two 
neighbouring sites. Two particular states of interest are 
the singlet 



FFa — 



FiF\ 



V2 



and triplet state formed from three states, 



FFa = 



F T F[ 



V2 



(14) 



(15) 



as well as F^ Fj and F± F[ . The two-site system described 
by all sixteen states of this form is not isolated. The 
rest of the system acts as a reservoir. The coupling is 
realized through the electron propagator for the rest of 
the system, which can be obtained from the appropriate 
trace of the single-particle Green function defined in Eq. 
(0) . Consequently, the energy levels of the two-site states 
are broadened. We define then the resolvent of a two- 
site state, R nm (t) = 9(t)((§ n (t)&l n (0))) , as a trace over 
the degrees of freedom of the reservoir and the two-site 
subsystem. The Fourier transform of the resolvent 



duj' 



oj — u' + iS 



(16) 



can be written in terms of the spectral function for the 
two-site system, which is defined as a mn = —lm.R nm /'K. 
The <5j's are then expanded 



Sj 



<J>t (f, 



(17) 



in terms of the sixteen level operators that describe the 
two-site physics. The products of Sj's that comprise the 
dynamical corrections are then expressed in terms of the 
two-site resolvents using the non-crossing approximation 
(NCA). Consequently, the total self-energy contains all 
the possible convolutions of the two-site resolvents. As 
the dynamics included are local, the NCA is expected 
to be accurate [Q. The resolvents associated with lo- 
cal triplet and singlet states are sharply peaked at well- 
defined energies Jl(J. The energy difference is small, how- 
ever, and given roughly by t 2 /U . Consequently, singlet- 
triplet mixing cannot be ignored. We can introduce this 
spin-fluctuation effect as a correction to the self-energies 
of the resolvents. The effective energy of the exchange 
interaction is given by the difference, 



J = 



Uj{(J FFs - (7 FFa )cL>, 



(18) 
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of the first moments of the spectral functions for the sin- 
glet and triplet two-site states. Antiferromagnetism cor- 
responds to J > 0. To facilitate a self-consistent eval- 
uation of the exchange energy, we computed the equa- 
tions of motion for the two-site level operators in the 
presence of the spin-fluctuation term, Jn • n Q /2. Here 
n = (ail , Cii)<Ji(cii , CiiY , with cr.; the Pauli matrices and 
a represents an average over the nearest neighbours of 
the second site in the cluster exclusive of the first site. 
A final quantity that is needed is the occupancy in the 
two-site states 

n$ n = Z$ n /Z, (19) 

where Z$ n = J d,Lua nn (uj), where Z represents a sum 
over all Z nm 's for the sixteen states and a nm — 
exp{-(3uj)<j nm (uj). 

III. RESULTS 

To obtain a self-consistent solution for the self-energy, 
we start with an initial guess for the electron spectral 
function which will be used to describe the properties of 
the environment in which the two-site system is placed. 
The two-site resolvents are then determined iteratively. 
Once the resolvents are determined, they are fed into the 
dynamical corrections and the Green function is deter- 
mined. A new spectral function is then computed thereby 
closing the self-consistent set of equations. This proce- 
dure is repeated until convergence is reached. 

A. Hubbard gap 

Shown in Fig. ([!]) is the total electron spectral function 
(-Im(Gn(fc, u) + 2Gi 2 (fc, u) + G 22 (fc, w))/tt) for the half- 
filled Hubbard model for U — 8t and T = 0.15*. Each 
trace corresponds to a momentum starting from (0, 0) to 
(-7T, 7r) to (tt, 0) and then back to (0, 0). Clearly shown are 
the upper and lower Hubbard bands with an energy gap 
of order U and the flatness of the band near the (tt, 0) 
point. The broadening is due entirely to the dynamical 
corrections. In the absence of such processes, the spec- 
tral function would correspond to a series of <5-functions 
at the lower and upper Hubbard bands. In the upper 
Hubbard band, the dominant spectral weight lies at the 
(tt, tt) point whereas in the lower band the spectral weight 
is located at (0,0). Broad spectral features obtain near 
the (tt, 0) point. The analogous spectral function for the 
D=l case is shown in Fig. (|J). As in the 2D case, the 
energy separation between the lower and upper Hubbard 
bands for D=l is striking. In both cases, the periodic- 
ity of the band is 2n as is expected for a paramagnetic 
solution. Nonetheless, we will show that local antiferro- 
magnetic correlations exist. 
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FIG. 1. Momentum and energy dependence of the electron 
spectral function for U = 8t and T — 0.15t for different val- 
ues of k. From top to bottom, the momenta correspond to 
(k x ,k y ) — (0,0) — > (7T,7r) — > (tt, 0) — > (0,0). Each momentum 
trace is shifted by hand. 

Integration of the spectral functions with respect to 
momentum yields the density of states (DOS) as a func- 
tion of energy shown in Figs. (||) and (Q). 
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FIG. 2. Momentum and energy dependence of the electron 
spectral function for U = 8t and T — 0.15t for different val- 



ues of k. From top to bottom, the momenta correspond to 
k x — — * tt. Each momentum trace is shifted by hand. 

It is evident that the Hubbard gap is fully formed for 
U fa W, where W is the bandwidth (W = At for D = 1 
and W = 8t for D = 2). In the D=l case, the gap is wider 
than in D = 2. For all values of U, we see a clear suppres- 
sion in the density of states near zero energy. However, 
for small values of U, the density of states at the Fermi 
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energy does not vanish at the temperatures we consider 
here. It is crucial then that we investigate the temper- 
ature dependence of the density of states at zero energy 
for small U. Shown in Fig. (||) is the density of states at 
the Fermi energy, p(0) as a function of temperature for 
both ID and 2D at U = 2t. 
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FIG. 3. Density of states for U = 8t,6t,tt,2t for T = .15* 
in D=l. The presence of a gap for all values of U indicates 
an absence of a metallic state at half-filling. 



For both D = 1 and D = 2, p(0) drops to zero as 
the temperature decreases. The enhancement seen in the 
density of states at zero energy in the D = 2 case over the 
D = 1 problem is tied to the shape of the non- interacting 
density of states. In D=l, as U decreases, the density of 
states becomes sharp (and actually diverges) at the band 
edges rather than approaching a constant value as for 
D = 2. In fact, in D = 2, the band exhibits a singularity 
at the gap edge. Consequently, it is easier to fill in the 
gap in D=2 than in D=l as is seen in Fig. (||). This 
trend also persists for any value of U as demonstrated in 
Fig. (§) for U = 8t. 
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FIG. 4. Density of states for U = 8t,6t,4t,2t for T = A5t 
in D — 2. The presence of a gap for all values of U indicates 
an absence of a metallic state at half-filling. 



We see clearly that p(0) is non-zero even in the strong- 
coupling regime, U ^> t, provided that the temperature is 
sufficiently large. This signifies that the physics at small 
and large U do not differ qualitatively. Because we have 
probed both U < W and U > W, we conclude that the 
Hubbard gap persists for all values of U both in D = 1 
and D — 2 and is certainly fully formed at T = 0. Hence, 
the only signular point in the half-filled Hubbard model 
is U = where the gap disappears. Consequently, there 
is an absence of a metallic state at half-filling. This result 
is in agreement with that of Moukouri and Jarell [jl2| and 
is consistent with the argument of Anderson || that the 
Hubbard model is always in the strong-coupling regime. 
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FIG. 5. Density of states at the Fermi energy for U — It 
as a function of temperature. 
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FIG. 6. Density of states at the Fermi energy for U = 8t 
as a function of temperature. 

Additional confirmation of the continuity between the 
weak and strong coupling regimes at half-filling is seen 
from the double occupancy shown in Fig. (Q). In terms 
of the two-site resolvents, the expression for the double 
occupancy is given by 

D = \ i n FDS + n FDA + n DB s + n D BA) + n DD (20) 
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where is given by Eq. (|l9|). In both cases, the dou- 
ble occupancy smoothly approaches the non-interacting 
value of 1/4 as U decreases. The absence of any kink 
in the double occupancy indicates an absence of a phase 
transition between the regimes, U < W and U > W. 
Hence, a continuity appears to exist between the small 
and large U regimes in the half-filled Hubbard models. 
The existence of a non-perturbative gap indicates that 
half- filled ID and 2D Hubbard models flow to strong 
coupling. 
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FIG. 7. Double occupancy as a function of U. In 
both D=l and D=2, the double occupancy approaches the 
non-interacting value of 1/4 as U — > and vanishes as 
U -> oo. 



B. Local Antiferromagnetism 

Because our method is capable of uncloaking the un- 
derlying spin dynamics, we investigated the behaviour of 
the local triplet, tiffs, and singlet, rip fa, occupancies 
for the states defined in Eqs. (15) and (|lj), respectively. 
These quantities were computed directly from the resol- 
vents that enter Eq. (|l^). Should the ground state be 
an antiferromagnet, the singlet occupancy, rip fa, should 
exceed the triplet occupancy, npFS at sufficiently low 
temperatures. For both ID and 2D the plots of the sin- 
glet and triplet occupancies are shown in Figs. (||) and 
(f|), respectively. 

As is evident, tiff a > tiffs as is consistent with an 
ordering tendency toward antiferromagnetism. Notice in 
the large U limit in 2D, singlet state formation is en- 
hanced over the ID value. However, to ascertain if these 
local probes are true signatures of ground state proper- 
ties, we computed the temperature dependence of tiffs 
and uffa for both ID and 2D at U = 8t. 




FIG. 8. Singlet (uffa) and triplet (tiffs) state occupan- 
cies as a function of U/t for D=l. For any non-zero value of 
U , uffa > tiffs indicating a tendency toward antiferromag- 
netic order. 




FIG. 9. Singlet (tiffa) and triplet (tiffs) state occupan- 
cies as a function of U/t for D—2. For any non-zero value of 
U, tiffa > tiffs indicating a tendency toward antiferromag- 
netic order. 



From Figs. ( |To| ) and (|TT|), we find thatat high temper- 
atures triplet excitations dominate. However, this trend 
is reversed below some temperature and the singlet occu- 
pancy becomes of order unity. This is significant and im- 
plies that the ground state is in fact an antiferromagnet. 
The tendency toward antiferromagnetism appears to be 
slightly enhanced in 2D relative to the ID problem. Our 
results then are consistent with antiferromagnetic order 
at T = in both ID and 2D. 
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FIG. 10. Singlet (jiffa) and triplet (uffs) state occupan- 
cies as a lunction of temperature for D=l. The fact that the 
singlet occupancy as T — > becomes of order unity is consis- 
tent with antiferromagnetic order. 
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FIG. 12. Effective exch ang e interaction J as a function of 
U/t computed using Eq. (tq). 



C. Heat Capacity 
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FIG. 11. Temperature dependence of the singlet (tiffa) 
and triplet {tiffs) occupancies for D=2. The fact that the 
singlet occupancy as T — » becomes of order unity is consis- 
tent with antiferromagnetic order. 



The energy splitting between the singlet and triplet 
states is due to an effective exchange interaction. Using 
Eq. (|l8|), we computed the effective exchange interaction 
shown in Fig. @ for both ID and 2D. Note first that 
J is always positive as a consequence of the fact that 
the singlet state is lower in energy than the triplet. This 
is a further indication of the antiferromagnetic order in 
the ground state. As expected, J is well approximated 
by t 2 /U in the strong-coupling regime. However, as U 
decreases, deviations from this behaviour are observed. 



The natural question that arises with any approximate 
treatment of strong correlation physics is: How seriously 
should the results be taken? Our study of the ID prob- 
lem is in part motivated by the fact that exact results are 
available from the Bethe ansatz. While Bethe ansatz is 
not amenable to yielding the Green functions from which 
the density of states can be calculated, ground state en- 
ergies and thermodynamic quantities are readily avail- 
able by this technique. Rather than compare with the 
total energy, we compute the temperature derivative or 
the heat capacity. Computation of the average energy 
is straightforward because we have already obtained the 
average double occupancy. Shown in Fig. (13) is a com- 
parison between the heat capacity computed within the 
present method (solid line) and the prediction from the 
Bethe ansatz |13 (triangles) for U = 8t. 
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1 

FIG. 13. Heat capacity for the ID half-filled Hubbard 
model as a function of temperature for U = 8t. The filled 
triangles are the results from the exact treatment via Bethe 
ansatz (see Ref. ( [D-3J)). 
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This figure demonstrates that at high to moderately 
low temperatures, the present method is quantitatively 
accurate, yielding results which differ by no more than 
1% from those of the Bethe ansatz. Such agreement is sig- 
nificant because in ID, correlation effects are particularly 
amplified. The two-peak structure of the heat capacity is 
tied to a competition between the contribution from the 
potential energy (high T) and the kinetic energy (low T) 
as illustrated in Fig. (O). 
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FIG. 14. Kinetic (ACt) and potential (ACu) energy con- 
tributions to the heat capacity of the ID and 2D half-filled 
Hubbard models for U — St. 



Both display maxima but in distinctly different en- 
ergy regimes. Near perfect agreement with the Bethe 
ansatz solution is obtained at high temperatures where 
the potential energy dominates. This is to be expected as 
the Hubbard operators provide an accurate treatment of 
the potential energy but only an approximate description 
of the kinetic energy. At sufficiently low temperatures, 
where the kinetic energy dominates, sharp spectral fea- 
tures appear and the numerical accuracy of the method 
wanes. Another source of error could be the two-site 
approximation itself. At low temperatures, an accurate 
description of the low-energy physics becomes essential. 
It might be that the two-site approximation inherently 
over-estimates the magnitude of the kinetic energy. To 
see if this breakdown persists for small U, we computed 
the heat capacity for U = 2t. The two-peak structure 
that occurs in the large U regime is absent for U -C W 
as illustrated in Fig. Jl5|). The disappearance of the 
two peaks is dictated by the non-interacting limit which 
possesses a single peak at T w Q.ht. Our results are in 
quantitative agreement with the numerical simulations 
of Shiba and Pincus |Q down to T sa At. Below this 
temperature, lack of numerical precision prohibited any 
accurate determination of the heat capacity. It appears 
then that the source of the breakdown at low temper- 
atures stems more from the lack of numerical accuracy 
than from the local description of the physics. However, 



more studies on this are necessary. 
0.7r 
0.6- 




° ' T/t 

FIG. 15. Heat capacity for the ID half-filled Hubbard 

model as a function of temperature for U = 2t. The devi- 
ation from the dashed line at low temperature stems from the 
emergence of sharp spectral features below T « Q.lt which 
prohibit an accurate numerical determination of the integral 
that enter the heat capacity. 



In 2D, the heat capacity (see Fig. (|l6|)) has the famil- 
iar 2-peak structure of the ID problem. Here again, this 
structure arises from a competition between the kinetic 
and potential energies as the dashed lines in Fig. ( |i"4| ) re- 
veal. Similar results have also been obtained by Scalettar 
and colleagues |2?j| from quantum Monte Carlo simula- 
tions on finite samples. The final feature on which we 
focus is the crossing of C(T, U) versus T for various val- 
ues of U . In ID comparison of Fig. ( jilt ) with (|l5| ) reveals 
that the heat capacities cross both at high temperature 
and at low temperature. Fig. ( [To] ) reveals that that this 
trend persists in 2D as well. The high temperature cross- 
ing point occurs at roughly T w 1.7t ± O.lt whereas the 
low-temperature crossing point is T ss At ± 0.1. The 
errors are due largely to the uncertainty in the data at 
small U. In the quantum Monte Carlo studies of Duffy 
and Moreo [Tj| on a 6 x 6 square lattice, similar values for 
the low and high temperature crossing points were found 
as well. In the D — > oo limit |l4|,[l(|, two crossing points 
are observed as well though at substantially smaller tem- 
peratures than in the 2D case. A unique crossing point 
for C(T, U) as a function of T for different values of U im- 
plies that at a particular temperature, the heat capacity 
is independent of U. This behaviour is observed in a wide 
variety of strongly-correlated experimental systems, such 
as 3 He JTJ, CeCu 6 - x Al x ft§, Nd 2 ~ x Ce x Cu0 4 |§, and 
UBe^ 1 20 1 . Vollhardt jl6| has shown that independence 
of C(T, U) on U at a particular temperature is fundamen- 
tally rooted in strong correlation physics. The condition 
for a unique crossing point for C(T, U) versus T for var- 
ious values of U can be recast ITfJ as 
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= 



dT <9C(T, U) 



T 



dU 



(21) 



in the limit that T — » oo. At high temperatures, 
C(T,U) oc U/T. Hence, dC/dU > as T -> oo. 
However, for the sum rule given by Eq. (|l|) to hold, 
dC/dU must change sign as the temperature is lowered. 
Such sign changes will be mediated by terms proportional 
to higher powers of U that enter with opposing signs. 
Hence, the sign change of dC/dU is a true correlation 
effect arising from terms at least proportional to U 2 and 
higher in the internal energy, E(T,U). As there is no 
phase transition as a function of temperature, the curves 
for C(T, U) must cross to satisfy the vanishing of the in- 
tegral in Eq. (pi]). At low T, the width of the crossing 
point is determined by low-lying excitations generated by 
the kinetic energy. The natural scale for such excitations 
is 4t 2 /U, in rough agreement with the low temperature 
crossing point in Fig. ([l6]). At high T, charge excitations 
dominate the contribution to the heat capacity. At large 
U, the gap should scale as U — W. 




T/t 



FIG. 16. Heat capacity for the 2D half-filled Hubbard 
model as a function of temperature for U = 2t, At, 8t. The 
crossing points at T m 1.7t ± O.lt and at T = At ± O.lt is in 
agreement with the general arguments of Ref. ( |l^|). The 
low-T peak in the heat capacity arises from spin fluctuations 
and the high temperature physics is tied to charge fluctua- 
tions. 



IV. CLOSING REMARKS 



We have applied a local method [ pL 0|Jl l| to the deter- 
mination of the dynamical corrections to the self-energy 
in the single-particle Green function for the ID and 2D 
Hubbard models. Although this method focuses on lo- 
cal 2-site correlations, it captures such established fea- 
tures as the onset of antiferromagnetism, absence of a 
Mott-Hubbard transition for non-zero U , and the uni- 
versal crossing of the heat capacity as a function of T for 
various values of U. As these features are the signatures 



of strong-correlation at half-filling, it appears that local 
dynamics offer an adequate description of these phenom- 
ena. Extending this method to 3 sites is prohibitive as 
this will entail an expansion in 4 3 = 64 three-site eigen- 
states. This calculation is impossible as the complexity 
of the two-site problem is already daunting. What does 
seem promising, however, is a possible field theory de- 
scription of the local dynamics that seem to be essential 
to an accurate description of strong correlation physics. 
Work along these lines as well as extending the present 
method to the doped case is underway. 
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